As a result of solving a system of linear equations, we obtain a vector of coefficients

\( u_{1,1} \\ u_{2,1} \\ u_{3,1} \\ \vdots \\ u_{k,l} \\ \vdots \\ u_{N_{x-2},N_y} \\ u_{N_{x-1},N_y} \\ u_{N_x,N_y} \\ \)

We remember that our solution (bitmap projection) is represented by a linear combination of two-dimensional B-spline functions

\( u(x,y) = \sum_{i=1}^{N_x} \sum_{j=1}^{N_y} u_{i,j} B^x_{i}(x) B^y_{j}(y) \)

In order to get the solution value in point \( (x,y) \) belonging to the area \( \Omega \) on which our bitmap is defined, we do not need, of course, to sum the values of all B-spline functions. We need to locate the element where the point \( (x,y) \) and all \( (p+1)^2 \) B-spline functions specified on this element. The other Bspline functions are zero at this point.

In our case, \( \Omega = [1,maxx]\times[1,maxy] \), and we have equally distributed \( N_x \) B-spline functions of order \( p \) along \( x \) axis and \( N_y \) B-spline functions of order \( p \) along \( y \). We have \( N_x-p \) elements along \( x \) and \( N_y-p \) elements along \( y \) axis. Each element has dimension \( [maxx / (N_x-p)]\times [maxy / (N_y-p)] \).

So to get the element number along the axis \( x \), we do division \( ne_x = int (\frac{x} {maxx/(N_x-p) }) \), and to calculate the element number along the y-axis, we do division analogously \( y \), we do division \( ne_z=int(\frac{y}{maxy/(N_y-p) } ) \), where int is the integer value (rounded down).

Then, on our element

\( [ne_x,ne_y] \) the functions are determined \( \{B^x_{i,p}(x)B^y_{j,p}\}_{ i=ne_x-p+1,ne_x+p-1,j=ne_y-p+1,ne_y-p-1 } \).

In order to calculate the solution value in point \( (x,y) \) we just need to calculate

\( u(x,y) = \sum_{ i=ne_x-p+1,ne_x+p-1,j=ne_y-p+1,ne_y-p-1 } u_{i,j} B^x_{i}(x) B^y_{j}(y) \)